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A re-investigation of the gravothermal catastrophe is presented. By means of a linear perturbation analysis, we study the dynamical 
stability of a spherical self-gravitating isothermal fluid of finite volume and find that the conditions for the onset of the gravothermal 
catastrophe, under different external conditions, coincide with those obtained from thermodynamical arguments. This suggests that 
the gravothermal catastrophe may reduce to Jeans instability, rediscovered in an inhomogeneous framework. We find normal modes 
and frequencies for the fluid system and show that instability develops on the dynamical time scale. We then discuss several related 
issues. In particular: (1) For perturbations at constant total energy and constant volume, we introduce a simple heuristic term in the 
energy budget to mimic the role of binaries. (2) We outline the analysis of the two-component case and show how linear perturbation 
analysis can be carried out also in this more complex context in a relatively straightforward way. (3) We compare the behavior of the 
fluid model with that of the collisionless sphere. In the collisionless case the instability seems to disappear, which is at variance with 
the linear Jeans stability analysis in the homogeneous case; we argue that a key ingredient to understand the difference (a spherical 
stellar system is expected to undergo the gravothermal catastrophe only in the presence of some collisionality, which suggests that the 
instability is dissipative and not dynamical) lies in the role of the detailed angular momentum in a collisionless system. 
Finally, we briefly comment on the meaning of the Boltzmann entropy and its applicability to the study of the dynamics of self- 
gravitating inhomogeneous gaseous systems. 
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1. Introduction 

Core collapse in an N-body system is a problem relevant to the 
dynamics of globular clusters, because these are considered to 
be the only stellar systems that possess the necessary degree of 
collisionality to relax thermally. The discovery of some clusters 
with cuspy cores was interpreted as a sign that they indeed expe- 
rienced the gravothermal catasttophe. In the context of globular 
clusters, the three main phenomena caused by collisionality are: 
core collapse, evaporation, and mass segregation. The focus of 
this paper is on core collapse. 

Consider the gravitational A^-body problem, where N » 1. 
That is, consider a set of N classical point masses, each of mass 
m, mutually interacting through Newtonian gravity. The particles 
may or may not be confined within a spherical volume of radius 
R. We are interested in the following questions: 

- Can the system reach some sort of equilibrium! Can this 
equilibrium be called thermal! 

- Can a model (for example, collisionless or fluid) of the N- 
body problem reach equilibrium? How does the equilibrium 
of a model relate to the equilibrium of the pure A^-body prob- 
lem and of a real stellar system? 

- Are these equilibria stable? 

This topic has been studied by many authors (for recent reviews 
see for example Heggie & Hut 2003 ; Bi nney & Tremainep 008). 
Different models of the Af-body problem (in particular gaseous 



models, collisionless models, and fluid models) admit equilib- 
rium configurations that are spatially ttuncated self-gravitating 
isothermal spheres, but it is not clear to what extent they can be 
considered as representative of the equilibrium states of the pure 
A^-body problem. Several studies have addressed the stability 
problem of isothermal spheres using thermodynamical methods, 
starting with |Antonov| ( |1962| > and |Lynden-Bell & Woodl(|1968|l 
and continuing with a long list of papers (for example, Hachisu 
| & Sugimoto| [19781 |Nakada| [T9781 |Katz||1978[ |Inagaki||198of 
Padmanabhan 1989; Chavanis 2002 2003 see also |Thirring| 
1 1970 1. In the thermodynamical approach the study of isothermal 
spheres is based on a form of entropy known as the Boltzmann 
entropy: 



S*[f] = ~k J"/(r,v)ln/(r,v)d 3 rd 3 v , 



(1) 



* Present address: Rudolf Peierls Centre for Theoretical Physics, 1 
Keble Road, Oxford OX1 3NP 



where / is the one-particle distribution function^jr and v the 
position and velocity vector respectively, and k the Boltzmann 
constant. 

Unfortunately, the thermodynamics of self-gravitating sys- 
tems still depends on a number of unresolved issues, partly be- 
cause of the long-range nature of the force and partly because 
of the divergent behavior of the force at short distance (e.g., see 
Padmanabhan 1990 ;lKatz|20031|Chavanis|2006t|Mukarrieil 2008; 
Cam pa et al.|2009||Bouchet et al.|2010 1. A critical analysis of the 
use of the Boltzmann entropy has been made by |Miller| ( |1973| l. 
We will briefly comment on this point in Sect. [2] 
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Therefore, it would be interesting to study the gravothermal 
catastrophe by limiting as much as possible the use of thermo- 
dynamical arguments. In this paper we reconsider the gravother- 
mal catastrophe and analyze the dynamical stability of self- 
gravitating fluids governed by the Euler equation by finding the 
normal modes and frequencies of spherical systems under vari- 
ous external conditions, in particular: constant total energy E and 
volume V (this case corresponds to the gravothermal catastrophe 
of Lynden-Bell & Wood 1968), constant temperature T and vol- 
ume V (isothermal collapse), constant temperature T and bound- 
ary pressure P (isobaric collapse, |Bonnor|[l956| Ebert |1955) . 
Note that, in contrast to thermodynamical arguments, the lin- 
ear modal analysis automatically gives the time scale for the de- 
velopment of the instability. The constant {T, V} case has been 
addressed by Semelin et al. ( |2001 1, who studied the stability of 
the system numerically, and by |Chavanis (2002), who found an 
analytical solution for the marginally stable perturbations us- 
ing methods developed by Padmanabhan (1989). The constant 
[T, P) ca se has been ad dres sed previousl y w ith synthetic argu - 
ments by |Bonnor| ( [T956l ) and |Ebert| | fl955) , by |Yabushita]([T 968), 
who studied the stability of the system numerically, by Lombardi 
& Bertin ( 2001) 1, who extended the analysis by Bon nor|(|1956| ) 
to the nonspherically symmetric case, and by |Chavanisf ( |2003[ ), 
who gave an analytical solution for the case of marginally sta- 
ble perturbations using the same method as he had used for the 
constant {T, V} case. The constant {E, V} case has been studied 
using a model based on the Smoluchowski-Poisson system (dif- 
ferent from the Euler-Poisson system considered in this paper), 
by |Chavanis et al.| ( |2002] ). 

In this paper, we extend these analyses of the constant {T, V} 
and {T, P} cases to the general calculation of eigenfrequencies 
and eigenfunctions for conditions outside those of marginal sta- 
bility, and study in detail the constant {E, V} case. We use a 
Eulerian or Lagrangian representation of hydrodynamics as sug- 
gested by the boundary conditions to be imposed on the system 
and we provide a unified treatment for all cases. Surprisingly, 
we find that the system becomes dynamically unstable in all the 
cases considered exactly at the same points found by |Lynden-| 
Bell & Wood ( 1968 ) by means of a thermodynamical analysis, 
that is for values of the density contrast (the ratio of the central 
density to the boundary density) as 14, 32. 1, and 709 respectively 
for the constant {T, P}, {T, V}, and {E, V} case. These results sug- 
gest that the gravothermal catastrophe may reduce to the |Jeans| 
(1902) instability, rediscovered in the inhomogeneous context. 
Then we briefly outline the problem of the linear dynamical sta- 
bility of a spherical fluid generalized to the two-component case, 
by perturbing the two-component spatially truncated isothermal 
sphere configurations previously considered by|Taff et al. ( 1975[), 
Lightman| ( |1977| ), |Yoshizawa et al.| ( |1978| ), |de Vega & Siebert| 
(2002 1, who studied the stability of these systems with a ther- 
modynamical approach, and Sopi k et al.| ( |2005| ), who addressed 
also the dynamical stability with a model different from the one 
used in this paper. We show that the onset of thermodynamical 
and dynamical stability occurs at the same point in the simplest 
case (constant {T, V}) and confirm that the component made of 
heavier particles is the primary driver of the instability. This re- 



sult is likely to be related to a similar finding by Breen & Heggie 
( |2012a|b| ) for two-component gravothermal oscillations. Finally 
we focus on the following puzzling phenomenon: in the colli- 
sionless model the instability disappears. The phenomenon is at 
variance with what happens in homogeneous systems (e.g., see 
Bertin 2000). In the last part of the paper we argue that the cause 
of the difference lies in the role of the angular momentum of the 
individual particles. 



The paper is organized as follows. In Sect.|2]we comment on 
the meaning of the Boltzmann entropy. In Sect. [3] we write the 
basic hydrodynamic equations in the Eulerian and Lagrangian 
representation. In Sect. |4] we show the results of the linear analy- 
sis of the hydrodynamic equations by studying the properties of 
spherically-symmetric perturbations. We find the relevant nor- 
mal modes and frequencies under different boundary conditions. 
In the analysis of the gravothermal catastrophe, we propose a 
modified expression for the energy, which can be considered as 
a simple way to incorporate the energy generation from bina- 
ries in the linear regime; we show that the catastrophe can in- 
deed be halted if we consider such a modified expression. In the 
last subsection and in Appendix [E] we briefly consider the two- 



component case. In Sect. 5. 1 we discuss the relevant time scales. 
In Sect. I5.2lwe discuss the difference between the collisionless 



model and the fluid model of the A^-body problem. In Sect. [6] 
we draw our conclusions and identify some open questions and 
issues. 



The thermodynamical approach and the 
Boltzmann entropy 



Antonov| ( |1962| ) and |Lynden-Bell & Wood] ( |1968| ) based then- 



stability analysis of self-gravitating spatially-truncated isother- 
mal spheres on the use of the Boltzmann entropy ([TJ. Starting 
from a kinetic description, they looked for stationary states of the 
Boltzmann entropy with respect to the distribution function /, at 
fixecQ total mass M, total energy E, and volume V = 4-nR 3 /3, 
where R is the radius of a spherical box. These stationary states 
are spatially truncated isothermal spheres. It is found that, for a 
given single-particle mass m, each spatially truncated isothermal 
sphere is identified by two dimensional scales (e.g., the central 
density po(0) and the temperature T) and one dimensionless pa- 
rameter {for example H = R/A, where A = [kT ImAnGpoiO)] ^ 2 is 
reminiscent of the Jeans length ( Jeans|1902 1; another equivalent 
choice for the dimensionless parameter is the density contrast 
po(0)/po(R), that is, the ratio of the central density to the density 
at the truncation radius R}. Conversely, each choice of the three 
quantities {T,p(0), E) identifies one particular isothermal sphere. 
This one-to-one correspondence may be lost for other choices of 
the quantities characterizing the system. 

In this thermodynamical approach the necessary condition 
for stability is that the stationary state corresponds to a (local at 
least) maximum of the Boltzmann entropy functional. |Antonov| 
( [1962} found that if S < 34.4 [i.e., if p (0)/p (R) < 709] the 
equilibrium configurations are local maxima, whereas for 3 > 
34.4 they are saddle points. Therefore, he concluded that self- 
gravitating isothermal spheres are unstable if S > 34.4. 

The work of Antonov ( 1962 1 was extended by |Lynden-Bell 
|& Wood| ( |1968[ ). These authors studied the thermodynamical sta- 
bility of the system under various conditions by applying the 
relevant thermodynamical potentials (based on the Boltzmann 
entropy) and gave a physical interpretation of the instability 
in terms of negative specific heats, as a characteristic feature 
of self-gravitating systems, thus creating the paradigm of the 
gravothermal catastrophe. They found that the critical value of 
S that corresponds to the onset of instability depends on the 



2 In terms of the distribution function /, we have M = m J/d 3 xd 3 v 
and 

„ f mv 2 s,3 r ,3 1„ 2 f/( r 'V)/(r',v') 3 3 3 , 3 , 

E= f(r,v)d rd v Gm drdvdrdv, 

J 2 JK 2 J |r-r'| 

where m is the single-particle mass. 
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adopted conditions. For example, S = 34.4 holds for a system at 
constant total energy E and volume V, S = 8.99 for a system at 
constant temperature T and volume V, and 3 = 6.45 for a system 
at constant temperature T and external pressure P. In this paper, 
we will show that the three points can be identified by means of 
a dynamical stability analysis of the fluid model. 

Now we briefly comment on the use of the Boltzmann ex- 
pression for the entropy ([1]). If we refer to systems for which 
the traditional thermodynamic limit (N — > oo and V — > oo at 
N/V fixed) is well defined and leads to homogeneous equilibria, 
it is well known ( |Jaynes|1965| l that the Boltzmann entropy cor- 
responds to the entropy defined in phenomenological thermody- 
namics only when the interparticle forces do not affect the ther- 
modynamic properties; that is, the Boltzmann entropy neglects 
the interparticle potential energy and the effect of the interpar- 
ticle forces on the pressure. If the equation of state is different 
from that of an ideal gas, the Boltzmann entropy is in error by 
a non-negligible amount. The correct expression for the entropy 
(corresponding to phenomenological thermodynamics) of sys- 
tems with a well-defined thermodynamic limit is the one given 
by Gibbs [for the definition of Gibbs entropy, see |Jaynes| ([T965 )]. 

The generalization of the above result to systems with inho- 
mogeneous equilibria, such as self-gravitating systems, is that 
the Boltzmann entropy is valid (i.e., it corresponds to entropy 
as defined in phenomenological thermodynamics) if and only if 
(1) the equation of state of the ideal gas applies locally to the 
equilibrium states; (2) hydrostatic equilibrium holds. Indeed, the 
most elementary way to understand the Boltzmann entropy in the 
case of self-gravitating systems is to think of a fluid in hydro- 
static equilibrium, with the equation of state of an ideal gas p = 
pkT/m, subject to reversible transformations between equilib- 
rium states (i.e., between different truncated isothermal spheres 
in the spherically symmetric case). As shown by Lynden-Bell & 
Wood ( 1968| > in their Appendix I, the classical thermodynamic 
entropy defined from the relation dS c = dQ/T and calculated 
from such transformations coincides with Sj,. For the general 
case of a system of particles interacting through an arbitrary 
two-body potential (not necessarily gravitational), it is possible 
to show that the stationary states of the Boltzmann entropy have 
the same density distribution as that of a fluid with the equation 
of state of an ideal gas in hydrostatic equilibrium. 

The above considerations suggest that the validity of the 
Boltzmann entropy is strictly related to the validity of the equa- 
tion of state of an ideal gas. In a kinetic description, the equation 
of state of an ideal gas is obtained by considering the local one- 
particle distribution function to be a Maxwellian (i. e., of the 
form f(r, v) — A exp[— mv 2 /2kT(r)], where A is a normalization 
constant) and by defining the pressure as p = J mf(v 2 /3) d?r d 3 v. 
This definition of pressure, which is obtained by considering par- 
ticles that would reverse their momentum when hitting an imag- 
inary wall, ignores the effects of interparticle forces. If, for ex- 
ample, particles repelling one another were confined in a box, 
they would exert some pressure on the walls of the box even if at 
rest: this contribution is completely neglected by the equation of 
state of the ideal gas (the neglected pressure is similar to that of 
rigid spheres when packed too closely; thus the Boltzmann en- 
tropy cannot give a correct result for a gas of almost rigid spheres 
when their density is too high). 

In the case of particles interacting through gravity, the ne- 
glected contribution is attractive and should therefore decrease 
the pressure compared to that of an ideal gas. Therefore, it is un- 
likely that the true thermal equilibrium state of an A^-body sys- 



tem when t — > oo (f is time) has an effective equation of state[jnot 
affected by the attractive nature of the gravitational force: thus 
the Boltzmann entropy may not be applicable to find the true 
thermal equilibrium of a self-gravitating A^-body system, since it 
assumes the equation of state of an ideal gas. Moreover, because 
gravitational forces have long range, each particle feels the in- 
fluence of all other distant particles. Thus, it may be that, strictly 
speaking, in the case of self-gravitating systems an equation of 
state cannot be defined in terms of local quantities. 

A different way to justify the use of the Boltzmann entropy 
is to show that it can be obtained from the Gibbs microcanoni- 
cal entropy by means of the so-called mean field approximation 
( |Katz|[2003j [Padmanabhan 1990). However, the range of appli- 
cability of this approximation is still not clear, also because a 
small-scale cut-off is necessary to avoid divergences in the mi- 
crocanonical entropy (Padmanabhan 1990). The Boltzmann en- 
tropy is also the H quantity of the Boltzmann //-theorem; from 
this point of view, a critical analysis of the Boltzmann entropy in 
the context of the gravothermal catastrophe has been performed 
by |MiUer| ( (T973) . 

The discussion in this section supports the hypothesis that 
self-gravitating isothermal spheres are not true thermal equi- 
librium states of the pure A^-body problem (i.e., the states that 
would be found in a spherical box containing N gravitating 
particles eventually, after an infinite amount of time), but are 
only metastable states of which the significance is still not 



completely understood (for example, see Padmanabhan 1990 
Chavanis 2006 and references therein). 



3. Basic equations of the dynamical approach 

In this section we derive the linearized hydrodynamic equations 
that govern the evolution of a fluid system for small deviations 
from the truncated isothermal sphere equilibrium configurations. 
We assume spherical symmetry, that is, we consider only ra- 
dial perturbations. The configurations are known to be stable 
against nonradial perturbations ( |Semelin et al.||200T) | Chavanis] 
|2l^|BTnney & Tremaine|2008) . 

Consider a self-gravitating fluid governed by the Navier- 
Stokes and continuity equations, together with the equation of 
state of an ideal gas (for example, see |Landau & Lifshitz|1987] l: 



^+V-(pu) = 0, 

at 



du V p 

— + (u • V)u = — - - VO + -V z u + 
at p p p 

kT 



P=P . 

m 

where p is the density, u is the fluid velocity, p is the local pres- 
sure, <t> is the gravitational potential, 77 and f are viscosity coeffi- 
cients, m is the one-particle mass, T is the temperature, t is time, 
and k is the Boltzmann constant. We will keep track of viscosity 
term until the end of Subsection |3.1| then (from Subsection |3. 2 
on), for simplicity, we will set 77 = £ = 0. In this paper, we do 
not perform an analysis of the effects of viscosity. In particular 
we do not discuss whether it has a stabilizing or destabilizing 
effect: hopefully, the role of viscosity will be studied in detail in 
a subsequent paper. 



3 By effective equation of state we mean an equation of state that, by 
imposing hydrostatic equilibrium, would reproduce the density distri- 
bution of the equilibrium state. 
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Under the assumption of spherical symmetry, the gravita- 
tional potential obeys the following equations: 



GM(f) 
V<5 = r--r, 



M(r) 



I 

Jo 



(3) 



p(s)4ns ds, 



which are equivalent to the Poisson equation. 

The hydrostatic equilibria of such fluid system are spa- 
tially truncated isothermal spheres (Chandrasekhar 1967) and 



are briefly described in Appendix [A As mentioned in Sect. [2] 
each self-gravitating truncated isothermal sphere is identified by 
two dimensional scales (for example the central density po(0) 
and the temperature T) and one dimensionless parameter S = 
R/A, which is the value of the dimensionless radius f = r/A 
at the truncation radius (the scale A was defined at the begin- 
ning of Sect. 2). Conversely, each choice of the three quantities 
{po(0), T, H) identifies a particular isothermal sphere. 

We shall denote unperturbed quantities by subscript and 
the perturbations by subscript 1. Thus the unperturbed density 
profile (the density profile of a truncated isothermal sphere) is, 
as a function of the dimensionless radius f (see Appendix [A): 



Po(f) 




= Po(0)e 



iff <E 
iff >a 



(4) 



where po(0) is the unperturbed central density and i^(f) is the 
regular solution to the Emden equation (the symbol ' denotes 
derivative with respect to the argument f): 



«A(0) = <A'(0) = 0. 



(5) 



When we perturb the hydrodynamic equations it is convenient 
to work in the Eulerian or Lagrangian representation of hydro- 
dynamics, depending on which boundary conditions we impose. 
In the Eulerian representation the independent variables are the 
time t and the position vector r. In the Lagrangian representation 
the independent variables are the time t and the original position 
ro that the fluid element under consideration had at the initial 
time fo. 

3.1. Eulerian representation 

In this subsection we derive the linearized perturbation equa- 
tions in the Eulerian representation. So we write each quantity 
as the sum of an unperturbed part (characterizing the truncated 
isothermal sphere) and a perturbation: 



p(r,f) =po(r) +pi(r,0 
u(r, t) — u\(r, t) 
T(t) = T + T x {t) , 



(6) 



where u is the radial component of the fluid velocity (recall that 
we consider only radial perturbations). We substitute Eq. |6]) in 
Eq. |2} and expand to first order in the perturbed quantities. We 
allow the temperature to vary in time while remaining uniform 
in space: this constraint will be used in order to impose the con- 
dition of constant total energy (see Subsection 4X2}. To find the 
normal modes of the system, we look for solutions of the follow- 
ing form: 

Pl {r,t) = p x {r)e- im 

u 1 (r,t) = u 1 (r)e- ,a " (7) 
7*1 (0 = fie~ m . 



In the following, we shall drop the symbol ' 



to keep the notation 
and us- 



simpler. After some manipulations (see Appendix B.l 
ing the unperturbed density profile Q we obtain the following 
linearized equation: 



Lf 



iu kTi _a 



4nGA m 
m I 1 



(8) 



10) 



Po(<Wo{ f 



Tlsf[?(fe*yY + (Z+i) 



where /(f) = po(f)Mi(f) is the unknown function, 



L = 



or 



4ttG Po (0) 



(9) 



represents the dimensionless (squared) eigenfrequency and £, is 
the following differential operator: 



£ = 



~df 2 



2 Ad 12 



T 



(10) 



Equation (|8]l is the equation to be solved to find the nor- 
mal modes and frequencies of the system, under the appropri- 
ate boundary conditions and constraints. The boundary condi- 
tions for Eq. ([8]) are defined in the following way. The absence 
of sinks and sources of mass, together with the assumption of 
spherical symmetry, implies /(0) = 0. Since in the Eulerian rep- 
resentation we shall consider only systems in a spherical box of 
constant volume, the radial velocity at the edge must be zero, 
which implies /(3) = 0. Thus the boundary conditions are: 



/(0) = /(H) = 0. 



(11) 



To solve the equations, we still need to specify the function 
T\(t). This specification discriminates between the constant en- 
ergy and the constant temperature case. Once this specification is 
made, Eq. ([H} is an eigenvalue problem: for fixed u>, the equation 
admits a solution only for discrete values of L. These solutions 
are the radial normal modes of the system. 

The operator £, has the following properties, which can be 
proved directly from the Emden equation Q: 



(12) 



These properties allow us to obtain analytical solutions in some 
cases. They are equivalent to those found by Padmanabhan 
(1989) in his review of the thermodynamical analysis of 
|Antonov| ( |1962[ ) and later also used by |Chavanis| ( |2002| |2003| l 
to obtain analytical solutions of the present hydrodynamic prob- 
lem fo r the constant {T, V) (Sect. |4. l.lj ) and constant {T, P) (Sect. 
|4.2.1| ) cases. 

Note that viscosity disappears when u> = 0, which is the sit- 
uation of marginal stability. Thus viscosity does not modify the 
points of the onset of instability (Sem elin et al.|200"Tj|Chavanis| 
120021). 



3.2. Lagrangian representation 

In this section we present the linearized perturbation equations in 
the Lagrangian representation in the inviscid case. In this repre- 
sentation, the independent variable is the position ro of the fluid 
element under consideration at the initial time fo- Thus in Eq. |2} 
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we have to perform the following change of independent van 
ables: 

[r -> r (r, t) 
\t -* t . 



(13) 



In the Lagrangian representation, each quantity is a function of 
the new independent variables ro and t. 



The calculations are summarized in Appendix B.2. 1 For lin- 
ear perturbations, the resulting continuity and Euler equations in 
the Lagrangian representation are (neglecting the term quadratic 
in the velocity): 



' d 


2 

r p 


d " 
dr 


P + 




at 


T*f>0 


r\ po dr \ 


' d 


2 

r p 


d N 

dr 0y 


it — 


r 2 dp 1 


dt 


r 2 p 


r\ dr p 



GM(ro) 



(14) 



As for the Eulerian case we separate each quantity in an un- 
perturbed part and a perturbation: 



p(r ,t) =po(ro)+pi(r ,t) 
u(r ,t) = ui(r ,t) 
r(r ,t) = r Q + n(r ,i) 



(15) 



We substitute Eq. ([15) in Eqs. ( p~4[ > and expand to first order in 
quantities with subscript 1. Then to find the normal modes we 
take: 

pi(r ,0 =pi('*o)e"^ 

Ul (r ,t) = Hi(r )e-''«" (16) 
n(r ,t) = ri(r )e- IW ' . 

In the following we shall drop the symbol " for simplicity of 
notation. After introducing the dimensionless radius £o = ro/A, 
using the density profile (|4) of the unperturbed state and after 
some manipulations we obtain the following equation for pi(£o) 
(see Appendix B.2.2 for an outline of the calculations): 



Pi 



-<!> d 



p„(0) ^ d^ 



jL e -*w r, 



± 



(17) 



where L = (J 1 /47rGpo(0) as for the Eulerian case. Boundary con- 
ditions are discussed in Subsection 4.2.1 Solving Eq. ( |T7] > al- 
lows us to find the normal modes in the Lagrangian representa- 
tion. We have allowed the temperature to vary in time in Eq. ( fP7) , 
but in the following we shall analyze only the isothermal case 
with Ti = 0. 



4. Modal stability of a self-gravitating inviscid fluid 
sphere under different boundary conditions 

Here we analyze the equations obtained in the previous section 
by imposing different boundary conditions. 



4.1. Eulerian representation 

In this subsection we analyze Eq. ([8]) by imposing two kinds of 
boundary conditions: constant temperature T and volume V or 
constant total energy E and volume V. 



4.1.1. Constant {T, V] case (isothermal collapse) 

Here we consider a fluid at constant temperature T contained in a 
sphere of fixed radius R. The condition of constant temperature 
is satisfied by imposing T\ — in Eq. ([8]). The condition of 
constant volume has been discussed in Sect. I3.ll and leads to 
the boundary condition /(E) = 0. Thus Eq. |8]l for the constant 
{T, V} case becomes: 

If = Lf (18) 

with the boundary conditions ( [TT| ). 

Equation ( fT8] l is an eigenvalue equation that, at given E, ad- 
mits solutions only for discrete values of L. If the lowest value 
of L at fixed S is positive, then all modes are stable (because u> 
is always real) and the system is stable. If the lowest value of L 
at a given value of E is negative, then unstable modes are present 
and the system is unstable. 

By means of the standard transformation 



/(0 = /(0exp 



(19) 



where £ is an arbitrary point in the domain of /, Eq. ( fl"8"j ) can be 
recast in the form of a Schrodinger equation: 



-f'(g) + f(g)U($ = Lf(&, 
where £/(£) is the effective potential given by: 

U <® = Ji + " f'^ ~ \ e ' m 

The boundary conditions become: 

/(0) = /(E) = . 



(20) 



(21) 



(22) 



The effective potential is shown in Fig. [T] From the boundary 
conditions (j22j> we see that choosing a specific S requires us to 
consider a potential that is infinite for £ > E. The existence of 
negative eigenvalues implies that the system is unstable. From 
the form of the potential it is clear that for small values of S neg- 
ative eigenvalues (i.e., unstable modes) do not exist. Negative 
eigenvalues appear only for sufficiently large values of S. For 
the present problem this occurs at S > 8.99. When S — > 00, an 
infinite number of negative eigenvalues appear, precisely at the 
same points where new unstable modes appear in the thermody- 
namical approach (see |Katz|1978| l. 

Let us now consider the solutions to Eq. ([18) in greater detail. 
Figure [2] shows the minimum value of L at given dimensionless 
radius E as a function of S. We see that L is negative, that is, 
the system is unstable, for E > 8.99, which is the same point 
found by |Lynden-Bell & Wood| (|T968) in the thermodynamical 
approach. Higher modes, that is, higher values of L at given S, 
would be represented by lines above the plotted curve. These 
lines would intersect the S axis at some points, which are the ze- 
ros of the analytical solution Gjy described below [see Eq. (|23)]. 
In Appendix D.l a few density and velocity profiles of numer- 
ically calculated eigenfunctions are shown. Most of them are 
for modes of minimum L at given E. Density profiles of higher 
modes exhibit oscillations not present in the lowest mode. 

For the case of marginal stability (L = 0), with the help of 
properties ( fT2) , the relevant eigenfunction can be expressed an- 



alytically (Chavanis 2002 



The function 



(23) 



is indeed a solution to Eq. ( [18) with L = 0, which satisfies 
Gtv(0) = 0. The values of S for which Grv(S) = are those for 
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Fig. 1. The effective potential U [see Eqs. ( |20| i and ( |2"Tj )] for the 
constant {T, V} case. To find the frequencies of normal modes, 
the Schrodinger equation |20} has to be solved with boundary 
conditions {22\ , which means that the effective potential is U (£) 
for £ < S and taken to be infinite for £ > H. States with nega- 
tive energies, which exist for S > 8.99, correspond to unstable 
modes. 



=8.99 



Fig. 2. The minimum value of the eigenvalue L at given S (the 
dimensionless radius characterizing the system) as a function of 
S, for the constant {T, V] case. If the minimum L is negative 
then the system is unstable. The system becomes unstable at 
S = 8.99, the same value obtained from the thermodynamical 
approach. From the figure we can read the typical time scale of 
the instability. As S — > oo, the asymptotic value L M = -0.042 is 
the same as for the other cases (see Figs. [3] and [4]!. 

which the boundary conditions ( fTT| are satisfied; thus they are 
the values of 3 at which each new unstable mode appears. The 
first zero of Gjv occurs where the first unstable mode appears, 
at S = 8.99. From the asymptotic behavior of iff it is possibile to 
obtain an asymptotic approximation of the zeros: they approxi- 
mately follow a geometric progression of ratio e 1 "^ ^ (see also 
|Semelin et al.|2001[|Cha7anis|20021 i. 

4.1.2. Constant {E, V} case (gravothermal catastrophe) 

Here we consider a self-gravitating fluid sphere at constant to- 
tal energy E and volume V. In this case the instability has 
been named gravothermal catastrophe by Lynden-Bell & Wood 
( 1968| l. The total energy E of the fluid is defined as: 



~NkT + „ 
2 2 



p(r)0(r)d 3 r. 



(24) 



It has two terms, which represent the thermal and gravitational 
contributions. The condition of constant energy is imposed in 
the following way. When the fluid is perturbed, its gravitational 
energy changes as a consequence of the redistribution of matter. 
We suppose that the temperature varies in time, while remain- 
ing uniform in space, so as to keep the total energy (thermal 
plus gravitational) constant (for a different model, based on the 
Smoluchowski-Poisson system of equations, the same nonstan- 
dard assumption has been made by Chavanis et al.||20 02). The 
thermal energy expression at time t is thus given by 3NkT(t)/2. 
In doing so, we are assuming infinite thermal conductivity (see 
also Sect. I5.ll for the relation of this fact to the relevant time 
scales). 

To reduce Eq. (|8]l to the constant {E, V} case we need to find 
the expression for the temperature as a function of the density 
distribution at fixed total energy, in the linear regime of small 
perturbations. Starting from Eq. ( |24| for the total energy and 
recalling that <I>(r) = -G JcFV p(r')/\r - r'|, we have: 



E = -NkT 
2 



_g r P (v)p{v' 

~2J |r-r'| 



} d 3 rdV. 



(25) 



By substituting T = 7q + T\ , p = po + p\ in Eq. ( |25| l, keeping 
only first-order quantities and imposing that the energy remains 
constant, we obtain the following expression for T\ : 



J Pl (r)(Do(r)d 3 r 



(26) 



Here (Do is the gravitational potential of the unperturbed density 
distribution po. After some manipulations (see Appendix \B.3\ 
we obtain: 



1 XW/(£)]/<W(£)d£ 
T i ~ ~— t — , r o 



io> §p (0MHV(H) 



(27) 



By substituting Eq. p7| ) in Eq. dSJ, recalling the definition A 
[kTo/4nGp[)(Q)m] 1 ^ 2 , and neglecting viscosity, we obtain: 



Lf = £f-tft 



— f 

SV(E) Jo 



<k#- [fm] d^ . (28) 



As discussed in Sect. 3.1 the boundary conditions are given by 



Eq. ( fTTj ). By integrating the last term of Eq. ( |28] i by parts under 
the boundary conditions ( fTT| , we obtain: 



Lf = Zf + ilf'e 



— f 

3 2 <A'(S) Jo 



(29) 



Note that Eq. p9|) [o r ((28J] contains an integral global con- 
straint. Equation (|29|l is to be solved to find the normal modes 
and corresponds to Eq. ( fT8| . The numerical procedure followed 
to solve Eq. ( [29) is relatively straightforward and thus is not re- 
ported here. 

Figure [3] shows the minimum value of L at given S. We see 
that the system is unstable for H > 34.36, which corresponds 
to a density contrast p (0)/p (E) > 709 ( |Antonov||1962| l. As H 
increases, new unstable modes appear, similarly to the behavior 
observed in the constant {T, V} case. As E — > oo, the asymptotic 
value of the minimum L is found numerically to be the same as 
in the constant {T, V} case, £«, -0.042. Such asymptotic value 
is the same also in the constant {T, P} case (see Se ct. |4.2. 1) >. 

In the case of marginal stability (L - 0), Eq. ( |29|l is equiv- 
alent to that found and solved analytically by Padmanabhan 
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1 a. 

'3 O 



==34.363 



-0.042 



Fig. 3. The minimum value of the eigenvalue L at given 3, the 
dimensionless radius characterizing the system, as a function of 
H, for the constant \E, V} case. When the minimum L is negative 
the system is unstable. The system becomes unstable at S = 
34.36, which corresponds to a density contrast =s 709, the same 
value obtained from the thermodynamical approach. From the 
figure we can read the typical time scale of the instability. As 
H — > oo, the asymptotic value of L, L x = -0.042, is the same as 
for the other cases (see Figs. [2] and [4). 



1989) in the thermodynamical approach. Here, for complete- 



ness, we record how the solution is derived, by adapting the 
method of Padmanabhan] (JT989J to our choice of variables and 
unknowns. 

To solve analytically Eq. ( |2"9"} for L — we rewrite it in the 
following form: 

Zf = -O'e^ , (30) 



where 



C 



= 1 f 

|E 2 (Zr'(E) Jo 



£V(0/(0df 



(31) 



is a constant (with respect to ^) that depends globally on /. Using 
the properties ( fT2} , we look for a solution of the form 



G EV (& = a^(0 + He 



(32) 



where a and b are real numbers. Since E q. (|30| > is linear in /, 
only the ratio a/b is relevant. Susbtituting (|32[> in |30} we obtain 
the first condition on a and b: 

a + b = C. (33) 

Using the expression of C ( |3T} and dividing by b, we obtain: 



a 
b 



i = 3 ~ 2 * MV® + € e ~*®) ^ ■ (34) 

|S 2 i/r'(S) Jo \b ) 



A second condition ona/b follows from the boundary condition 
/(H) = 0: 

V(S) + Ee-^ E > = . (35) 
b 

Substituting a/b from ( |3"5j ) in ( [34| > we obtain: 
3 



(36) 

The values of S satisfying Eq. ([36} are those for which 
Eq. ( [29} admits a solution for L = 0. In particular, the lowest 
value of 3 which satisfies Eq. ((36} determines the threshold of 



instability: by numerically solving the algebraic equation ( |36} , 
this minimum value is found to be S = 34.36 (as shown by 
Antonov 1962). As for the constant {T, V} case, the number of 
unstable modes for each S coincides with the results of the ther- 



modynamical analysis [see |Katz (1978}]. Once a value of S is 
obtained, the value of a/b and then an analytical solution Gev is 
determined. 

Density profiles for modes of minimum L at given H are 
shown in Appendix D.2 Since the equations involved are equiv- 
alent, the density profile for the marginally stable perturbation 
(L = 0) is the same as that found by Padmanabhan ( 19891. He 
pointed out that it has a "core-halo" structure, that is, an oscil- 
lation in p\ \ the density perturbation is positive in the inner part 
(nucleus), negative in the middle (emptying area), and then posi- 
tive again (halo). The core-halo structure has been physically in- 
terpreted in the framework of the gravothermal catastrophe given 
by |Lynden-Bell & Wood| ( |1968[ ), using the concept of negative 
specific heat. However, this interpretation is not applicable to the 
present context. 

We found that for modes of minimum L at given 3 the core- 
halo structure is present if L < 0.021, disappears between L - 
0.021 and L = 0.022, and is absent for L > 0.022, as illustrated 



in the figures presented in Appendix D.2 Higher modes always 



exhibit one ore more oscillations, as in the constant {T, V} case. 

We have thus shown that for the present case the behavior of 
the eigenvalues Las a function of 3 is similar to that of the con- 
stant {T, V} case, with the difference that the instability threshold 
is higher in the constant {E, V} case. The interpretation of this 
fact in the context of the fluid model is as follows. When the fluid 
is compressed, the gravitational potential energy decreases and 
the temperature increases in order to maintain the total energy 
constant. Therefore, the tendency toward collapse is weakened 
and instability can take place only at higher values of 3 (with 
respect to the case in which the temperature remains constant). 

In fact, if instead of expression ( [24} for the total energy we 
consider a modified expression in such a way that the temper- 
ature increase is greater, the collapse can be halted completely. 
Consider the following heuristic expression for the total energy: 



~NkT + „ 
2 2 



p(r)<D(r)dV 



vo- 2 R 3 p(0), 



(37) 



which contains an additional term proportional to the dimension- 
less parameter v and to the central density p(0). The quantity 
<Tq = kTo/m is the unperturbed thermal speed. We repeated the 
analysis of this subsection with such a modified expression for 
the energy and found the new value of 3 for the threshold of the 
instability. In this analysis the expression for T t ([27} obtained 
previously, to be substituted in Eq. (|8}, is replaced with the ex- 
pression for T\ that is obtained from Eq. ( [37} by substituting 
T - Tq + T\,p = po+Pi an d by imposing that the variation of the 
total energy E is zero. We found that when v > (i.e., when the 
fluid is compressed the new term contributes to make the tem- 
perature increase) the (linear) instability is postponed to higher 
values of 3 for small v and is completely halted for v > 5 x 10~ 4 . 
We also verified that if v < (i.e., when the fluid is compressed 
the new term contributes in the opposite direction), instability 
occurs at lower values of 3. 

In the case of the pure A^-body problem, it is believed that 
collapse can be halted by energy "generation" through bina- 
ries, giving rise to a phenomenon called gravothermal oscil- 
lations (for a review, see |Heggie &"Tl ut 2003 ). Gravothermal 
oscillations were discovered by Sugimoto & Bettwieser (1983) 
using a gaseous model. These authors introduced in the model 
of |Lynden-Bell & Eggleton (1980) a phenomenological energy 
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generation term (to represent the role of binaries), proportional 
to a power of the density and found that collapse can be halted 
and reversed. Similarly, our v term shows in a simple manner 
that the instability can be halted in the linear regime with a suit- 
able term in the total energy budget that mimics effects presumed 
to be associated with the presence of binaries. 

Gravothermal oscillations have been confirmed by iV-body 
simulations (Makino 1996). We will comment further on this 
topic in Sect. 5.2" 



<5 0.2 ■ 



4.1 .3. Constant {T, V}: two-component case 

We briefly studied the problem of dynamical stability by means 
of a linear modal analysis of a two-component ideal fluid. Each 
component is assumed to interact with the other only through 
the common gravitational potential. The hydrostatic equilibrium 
configurations are spatially truncated two-component isothermal 
spheres, as considere d by |Taff et al. ( 1975 1, Lightman1(|1977[l, 
and Yoshizawa et al. ( 1978 1; see also de Vega & Siebert] ( |2002| i 
and |Sopiketal.| ( |2005| l. 



We call the single-particle masses ma and mj, with mj > m^. 
Two-component isothermal spheres are characterized by two 
additional dimensionless parameters (with respect to the one- 
component case): the ratio of the single-particle masses m^/m^ 
and the ratio M B /M A of the total masses associated with the two 
components. The reader is referred to Appendix [E] for a descrip- 
tion of the equations used. 

We considered the constant {T, V} case, which is the sim- 
plest from the mathematical point of view, and found that the 
equivalence of the thermodynamical and dynamical approaches 
still holds. It is possible to show analytically that the onset of dy- 
namical instability takes place at exactly the same points as those 
found with the thermodynamical approach. The analysis in the 
thermodynamical approach was performed by generalizing in a 
straightforward manner the analysis of Chavanis (2002). 

An interesting result of this analysis is that the instability ap- 
pears to be driven by the heavier component, in the following 
sense. Consider the ratios PxaIpq and pislpa of the density per- 
turbation of each component to the total local unperturbed den- 
sity. The ratio referring to the heavier component can be higher 
even if the total mass of the heavier component is very small. 
For example, for m B /m A = 3, we found that the ratios PiaIp®, 
Pib/po were comparable for Mb/Ma - 0.066 (see Fig. E.li. 
For higher values of Mb/Ma the heavier component dominates. 
Independent indications that the heavier component dominates 
the collapse have been found by Sopik et al. (2005 i for a differ- 
ent model based on the Smoluchowski-Poisson system of equa- 
tions. Breen & Heggie (2012a b) found that the heavier com- 
ponent dominates gravothermal oscillations in two-component 
clusters. This is likely to be related to the present analysis. 



4.2. Lagrangian representation 

4.2.1. Constant {T,P} case (isobaric collapse) 

In this section we analyze perturbations at constant temperature 
T and constant boundary pressure P. Therefore, T\ = and 
Eq. ( fTT) becomes: 



Pi 
Po(0) 



-<A(fo) J 



«0 



(38) 



Fig. 4. The minimum value of the eigenvalue L at given H, the 
dimensionless radius characterizing the system, as a function of 
S, for the constant {E, V} case. When the minimum L is neg- 
ative the system is unstable. The system becomes unstable at 
S = 6.45, the same value as found in the thermodynamical ap- 
proach. From the figure we can read the typical time scale of the 
instability. As H — > oo, the asymptotic value L x = -0.042 is the 
same as for the other cases (see Figs. [2] and [3J. 



Boundary conditions are as follows. The condition mi(0) = 
translates into the condition p\(0) - [using Euler equation 
( |14) , under the assumption that du/dro does not diverge and 
du/dt = -icuu]. The condition of constant pressure requires that 
a fixed Lagrangian fluid shell (which follows the fluid during the 
motion and thus does not have a fixed position in space) feels 
constant pressure. Constant pressure requires constant density 
because of the ideal gas equation of state. Therefore, the rele- 
vant boundary conditions are: 



Pi(°) : 
Pi(H) 





:0. 



(39) 



From mass continuity and by imposing u\ (0) + (which is true 
in the Eulerian representation and thus we expect to be true in 
the present case), we also have the condition p\(Q) + 0. 

The numerical procedure to solve Eq. ( |38") > is relatively 
straightforward and is thus not reported here. Figure ^repre- 
sents the minimum value of L at given E. The system becomes 
unstable for E > 6.45, which is the same condition as found by 
Bonnor| ( [T936) , [Ebertl ( [T955) , and |Lynden-Bell & Wood| ( p68) . 
As H increases, other unstable modes appear, similarly to the 
cases described previously. As S — > oo, the asymptotic value of 
the minimum L is found numerically to be the same as for the 
constant [E, V} and {T, V} cases, L M = -0.042, but, in contrast 
to the constant {E, V} and {T, V} cases, it is reached from below. 
There is a minimum at L w -0.043 for H ~ 15. 

In Appendix |D.3| density profiles of normal modes are 
shown. Most of them are for modes of minimum L at given 
S. Higher modes present oscillations. It would be interest- 
ing to show whether the analysis of this subsection might be 
generalized to the nonspherically-symmetric case by following 
|Lombardi & Bertin | ( |200T] ). 

As for previous cases, we can obtain analytical solutions for 
the marginally stable perturbations. For L — 0, Eq. (|38]l reads: 



Pi 
Po(0) 



(40) 
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where X is defined as: 



£ = -! + 



e~* d 



go 



4_ 



(41) 



From Emden Eq. Q, the operator X. has the following proper 
ties: 

£0A' 2 ) = -^ 



(42) 



Hence, an analytical solution of Eq. ( |40[ > is the following (for 
a derivation in the Eulerian representation, see also Chavanis 
[2003) : 

G T p(^) = 2e~* - <A' 2 . (43) 

This solution satisfies the boundary conditions pi(0) = 2 and 
p'j(0) = 0, as required by Eq. ([39). The zeros of Gtp(^q) allow 
us to identify the marginally stable normal modes that satisfy the 
correct boundary conditions ( |39) . The first zero of Gtp(^o) is at 
S = 6.45. Other zeros correspond to the values of S at which 
new unstable modes appear and are found to be the same as in 
the standard thermodynamic al approach. 

5. The different behavior of a collisionless 
self-gravitating sphere 

5.1. Time scales 

In this subsection we briefly discuss the typical time scales that 
characterize gaseous and fluid models and compare them to 
those of globular clusters. To a large extent, our discussion fol- 
lows that of Inagaki]( [T980) . 

Let Tlte be the local relaxation time, tqje the global relax- 
ation time, and tj the dynamical time scale. For a gaseous sys- 
tem, we identify Tlte with the time (of the order of the inverse 
mean collision frequency) needed to reach a local Maxwellian 
distribution, t cte with the time in which thermal conduction 
balances the temperatures of different parts of the system, and 
Td with the sound travel time. Gaseous models generally assume 
the following ordering: 



T LTE ^ T d ^ T GTE 



(44) 



For a globular cluster, because the mean free path is very large 
and stars can cross the system many times before actually collid- 
ing, the ordering of time scales is different. Stars do not collide 
significantly with neighboring stars, following the mechanisms 
that usually characterize thermal conduction; instead, they tend 
to release their energy through the entire cluster (also because of 
the long range nature of the force). Thus for globular clusters we 
have: 

Td « T LTE =i t cte . (45) 

In the fluid model analyzed in this paper we assumed infinite 
thermal conductivity, because the temperature was always taken 
to be and to remain uniform. Moreover, we implicitly assumed 
that the distribution function is locally Maxwellian. Therefore, 
our fluid model follows the ordering: 



T LTE Td 
T GTE Td ■ 



(46) 



Note that assumptions ( |44) an d (|4"6| ) are not applicable to the 
situation of a globular cluster (|45|l. This difference and the re- 
lated limitations must be kept in mind if we wish to apply these 
models to understand the evolution of globular clusters. 



5.2. The dynamics of a collisionless self-gravitating 
isothermal sphere 

Isothermal spheres are equilibrium configurations for several 
idealized models of the pure A^-body problem. In particular, self- 
gravitating isothermal spheres can be studied as stationary states 
of the collisionless Boltzmann equation: 

a/(r,v,Q +v a/(r,v,Q dO(r,Q df(r,v,t) =Q 



dt 



By 



By 



Therefore, it is natural to ask whether such collisionless isother- 
mal spheres are or can be unstable. It can be shown ( |Binney| 
& Tremaine 2008 ) that the unbounded collisionless isothermal 
sphere is linearly stable, in contrast to the results of the fluid 
counterpart (in the fluid model we recover the unbounded case 
by taking H — > oo). It is generally believed that the collisionless 
isothermal sphere is also stable in the nonlinear regime, although 
to our knowledge a rigorous proof of this statement is still lack- 
ing. Moreover, if only spherically symmetric perturbations are 
considered, it is possible to show, through an argument based on 
the conservation of the detailed angular momentum^] (Appendix 
|C) , that a collisionless sphere (bounded or unbounded) cannot 
collapse, because each particle has a minimum radius it can at- 
tain. In contrast, the fluid system analyzed in Sects. [3]and[4]is un- 
stable only with respect to spherically symmetric perturbations, 
with the subsequent nonlinear evolution presumably leading to 
a collapse. The above considerations support the hypothesis that 
self-gravitating collisionless isothermal spheres (bounded or un- 
bounded) are stable with respect to all kinds of perturbations and 
cannot coll apse ( |Kandrup & Sygnet||1985[ |Kandrup||1990[ |Batt 
|etal.|1995| >. 

The different behavior, between the collisionless and the 
fluid case, emerged in this paper is actually quite surprising, be- 
cause for the homogeneous case the collisionless and fluid mod- 
els behave in the same way with respect to Jeans instability (e.g., 
see Bertin 2000), in the sense that both the fluid and the colli- 



sionless systems are linearly unstable under the same criterion 
for instability. 

If the self-gravitating collisionless isothermal sphere were 
unstable, its instability would develop on the dynamical time 
scale; in turn, it is commonly believed that the gravothermal 
catastrophe may occur only if the system is at least weakly col- 
lisional and thus it is thought to develop on the collision time 
scale. 

In the spherically symmetric case a difference between the 
two models is the following: while in the fluid model each fluid 
element is sustained against gravity by pressure of the inner 
parts, in the collisionless model stars are sustained by their indi- 
vidual angular momentum relative to the center, that is, by their 
velocity dispersion]^] In moving from a kinetic description to a 
fluid description, all the information about microscopic veloci- 
ties and detailed angular momentum is lost: in particular, each 
fluid element has zero angular momentum (see also Appendix 
O. In this respect, the real situation of a globular cluster resem- 
bles more the collisionless case: strictly speaking, stars are not 
sustained by pressure, but rather by their velocity dispersion, be- 
cause the mean free paths are long and stars cross the cluster 
many times before feeling the effects of collisions. 



4 By detailed angular momentum we mean the angular momentum of 
the individual particles. 

5 Even if the total angular momentum vanishes, the sum of the mag- 
nitudes of the angular momenta of the individual particles is different 
from zero. 
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If we come back to the original pure A^-body problem, it is 
therefore natural to ask: what is in this case the role of detailed 
angular momentum? We may argue that a quantity related to the 
detailed angular momentum should characterize the process of 
core collapse. If the collapse can happen in an A^-body system, 
it may be accompanied by a slow decrease of the sum of the 
magnitudes of the angular momenta of the individual particles 
slowly sinking toward the center. 

For a particle with angular momentum / that moves in a 
gravitational potential generated by a mass M a radius related 
to angular momentum (which is the radius of the circular orbit if 
M is a point mass) can be defined as: 



d = 



J 1 



Gm 2 M ' 



(48) 



Therefore, for a stellar system of N stars and total mass M we 
can define a radius related to detailed angular momentum in the 
following way: 

A 



^ am ~ NGm 2 M' 



(49) 



where A is the sum of the squares of the angular momenta of the 
individual stars, that is, the following quantity: 



(50) 



Thus we may argue that the typical time scale for core col- 
lapse should correlate with: 



tec — 



dA/df 



(51) 



The main mechanism through which A varies with time is ex- 
pected to be that of two-body collisions. Thus f cc should be of 
the order of the two-body relaxation time. 

By monitoring the quantity A(f) in A^-body simulations, it 
would be interesting to test the role of detailed angular mo- 
mentum in the mechanism of gravothermal oscillations (Makino 
1996| > [A(?) may reach an equilibrium value, around which the 



system gravothermally oscillates; this would be consistent with 
the fact that the typical time scale of gravothermal oscillations 
is the two-body relaxation time], its relevance to the studies of 
core-collapse in the gaseous model ( |Lynden-Bell & Eggleton 



1980| |Sugimoto & Bettwieser|[1983[ ), and its connection with 
the phenomenon of the gyro- gravothermal catastrophe ( |Hachisu| 



6. Discussion and conclusions 

In this paper we have studied systematically the dynamical sta- 
bility of a self-gravitating isothermal fluid sphere by means of a 
linear modal analysis with respect to spherically symmetric per- 
turbations. In this sense, we have studied the Jeans instability 
in the inhomogeneous context of a sphere of finite size. Within 
a unified framework, by imposing the boundary conditions of 
constant {T, V) (isothermal collapse), constant [E, V) (gravother- 
mal catastrophe), and constant [T, P) (isobaric collapse), we 
have proved that the onset of dynamical instability occurs ex- 
actly at the points identified in the thermodynamical approach 
(see|Antonov|[T962l|j^nden-Bell & Wood|1968||Bonnor|1956| 
Ebert ||1955 1 and, by adapting derivations from other authors 
Vadmanabhanl[T990l |Chavanis||20"02] [2003) , we have provided 
an analytic expression for the eigenfunctions of the marginally 



stable modes. Indeed, as noted in the Introduction, some results 
along these lines have been obtained previously by other authors. 
The main new results obtained in this paper are the following: 

- Using the fluid model based on the Euler equation, we have 
extended previous studies of the constant {T, V} and {T, P} 
cases to the constant {E, V} case (gravothermal catastrophe), 
proving that the onset of Jeans instability occurs exactly at 
the same point identified in the thermodynamical approach 
also in this case. For this constant {E, V} case, we have in- 
troduced a heuristic term to incorporate effects akin to the 
stabilizing role of binaries. 

- For all the three cases described above, we have calculated 
numerically eigenfrequencies and eigenfunctions of the rele- 
vant modes also outside the conditions of marginal stability. 
The time scale for the instability that we have found is the 
dynamical time scale. The excitation of higher modes has 
been illustrated in a simple way, by referring to an effective 
potential that governs the structure of the linear modal anal- 
ysis. 

- We have found that for all the cases treated in our investi- 
gation, as the dimensionless radius of the isothermal sphere 
H becomes larger and larger, the value of the dimensionless 
growth rate of the most unstable mode tends to a univer- 
sal asymptotic constant value, independent of the adopted 
boundary conditions. 

- We have briefly shown that the correspondence between 
the stability in the dynamical and in the thermodynamical 
approach also holds for the two-component case and have 
found indications that the heavier component is the more im- 
portant driver of the instability. 

- As a general discussion, we have commented on the mean- 
ing and applicability of the Boltzmann entropy for self- 
gravitating systems and argued that the main difference be- 
tween the dynamical behavior of a fluid and a collisionless 
sphere, in relation to their application as models of the pure 
A^-body problem or a real weakly collisional stellar system 
such as globular clusters, is to be ascribed to the role of the 
detailed angular momentum behavior in the collisionless and 
weakly collisional cases. 

The role of the viscosity and its consequences outside the con- 
dition of marginal stability have not been examined; hopefully, 
this issue will be addressed in a future paper. 

Another interesting question is how the instability depends 
on the particle-particle interaction, for non-Newtonian cases 
(Padmanabhan 1989). Potentials that exhibit a softe ning at small 
radii, such as 1 / (r L + r^) 1 ^ 2 , where ro is a constant ( [Chavanis &| 
Ispolatov 2002 see also Casetti & Nardini|2012 1, or that decline 
with a different power law at large radii, such as 1 / r", with a + 1 
( jlspolatov & C ohen 2001), have indeed been considered. Based 
on the present article, we may argue that the results obtained 
from the thermodynamical approach would be reinterpreted and 
clarified as the analogue of the Jeans instability, by studying a 
fluid model for the potential considered, with the equation of 
state of a perfect gas. 

In general, this paper strengthens the view that the applica- 
bility of different idealized models to describe the process of 
core collapse in systems made of a finite number of stars is more 
subtle than commonly reported and that, in general, the study of 
Jeans instability of inhomogeneous stellar systems still leaves a 
number of questions open. 

Acknowledgements. We would like to thank Marco Lombardi, Francesco 
Pegoraro and Steven N. Shore for many interesting comments and discussions. 
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poW = pa(f)kT/m 



(A3) 



We express M(r) by means of the condition of hydro- 
static equilibrium ( [A.l[ ), differentiate to obtain dM, and equate 
it to dM = 4npo(r)r 2 . By making the change of variable 
Po(r) = po(Q)e^ (r \ where po(0) is the central density, and 
introducing the dimensionless radius £ = r/A, where A = 
[kT /47rGpo(0)'w]^ 2 , we obtain the differential equation for ifr(^), 
recorded in the main text as Eq. d3). [Because the constant 
po(0) is interpreted as the central density, we are considering 
the boundary condition i//(0) = 0; the density is taken to be 
regular at the origin, so that the second boundary condition is 
tfr'(Q) = 0.] The solution of Eq. |5]l (called Emden equation) is a 
monotonic increasing function characterized by logarithmic be- 
havior ifr(g) ~ ln£ 2 and iA'(^) ~ 2/£ as £ — > oo . 

From Eq. |5]) the mass enclosed within the radius £ is: 



kTA 
Gm 



£V(£) 



(A.4) 



From Eq. ( |A.4| i and the asymptotic behavior of ip, it is clear that a 
solution with finite total mass is obtained only by truncating the 
system at a dimensionless radius £ = 3. A truncated isothermal 
sphere is then identified by two scales T and po(0) and one di- 
mensionless parameter S. The density profile of a spatially trun- 
cated isothermal sphere is given by po(r) = po(0)e~^ r) where i/r 
is the solution to Eq. d5). 

Appendix B: Linearization of the hydrodynamic 
equations 

6.7. Eulehan representation 

Here we record the calculations leading to the linearized Eq. ([8). 
The unperturbed density profile is given by Eq. Q. We substitute 
Eqs. (|6]l in Eq. (|2} and expand to first order in quantities with 
subscript 1 to obtain: 



dpi 
dt 



+ V-(p Ui) = : 



dui _ / Vp pj Vp, \ kT Vp () kT 
dt 



Po Po Po / m 



■ AnG 



(B.l) 

J^pi(s,t)s 2 ds 



Po m 



C + - 

" V 2 ui + -V(V-m) 



Po 



Po 



(B.2) 

Then we look for solutions of the form |7]). From Eq. ( |B.1| > we 
obtain p\. 

V-(p Ui) 

Pi = : , (B.3) 



and thus eliminate it from Eq. ( |B.2| i to find the radial component 
of the Navier-Stokes equation: 



Appendix A: Fluid truncated isothermal spheres 

In this appendix we summarize the properties of spatially trun- 
cated self-gravitating fluid isothermal spheres. 

The equations for the hydrostatic equilibrium of a spherically 
symmetric fluid with the equation of state of an ideal gas are: 



GM(r)p (r) dp Q (r) 



M(r) 



f 

Jo 



dr 



4ns p Q (s)ds . 



(A.l) 



(A.2) 



ldpo/dr 1 d 

u) p ui = 



(r 2 po«i)-^- 4^-(>- 2 P(>mi) 
po r L or v ' or [r z or y ' 



kTo 
m 



. dp Q kT 1 2 

-iu>— 47rGp (1 Mi 

or m 

.[Id/ 2 dui\ r] d 



r 2 dr \ dr 



3 dr 



\_a_ 

r 2 dr 



ir 2 m) 



(B.4) 



where U\ is the radial component of the velocity. By defining 
f(r) = po(r)u\(r) and introducing the dimensionless radius ^ = 
r/A we obtain Eq. ([8]). 
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B.2. Lagrangian representation 
B.2.1 . Change of variables 

Here we show the change of variables leading from Eqs. Q 
(Eulerian representation) to Eqs. ([14} (Lagrangian representa- 
tion). Let us assume spherical symmetry and neglect viscosity. 
By dropping the nonlinear term (u ■ V)u, the Euler equation and 
the continuity equation |2} become, in the Eulerian representa- 
tion: 

d P _l_ 1 d t 2 \ n 
ft + 7*dr (r PM) = °' 

du _ dp/drkT GM(r) (B ' 5) 
dt p m r 2 

Now we perform a change of variables. In the Lagrangian rep- 
resentation, each quantity is described as a function of the new 
independent variables ro and t, where ro is the position of the 
fluid element at t — to- The standard rules to transform deriva- 
tives of a generic function / are: 



df(r,t) dr (r,t)df(r ,t) 



dr 
df(r,t) 
dt 



dr 



dr 



dr (r,t)df(r ,t) df(r ,t) 



(B.6) 



dt 



dr 



8t 



The partial derivatives of ro are obtained from the following re- 
lation, which expresses the condition that two fluid shells do not 
cross each other: 



Jrro(r,t) 
PO 




po(s)47rs 2 ds = I p(s, t)4ns 2 ds 



(B.7) 



By taking the partial derivative with respect to r of Eq. ( |B.7| > 
(each side of the equation is considered as a function of r, f) we 
obtain: 

2 , .dr (r,t) 2 



r po(ro)- 



dr 



= rp{r,t) 



(B.8) 



By taking the partial derivative of Eq. ( |B.7| i with respect to t we 
obtain: 

r Po(ro) ^ — = I * — ET~ ds ■ ( B9 ) 



dt J,, dt 
From the continuity equation, Eq. (|B.9|l then becomes: 



roPo(ro) dr °g t '^ = -r 2 p(r,t)u(r,t) 



(B.10) 



From Eqs. ( |B.10| i, ( |B.8| l, and ( |B.6| l we obtain the equations of 
hydrodynamics in the Lagrangian representation ( fT4] >. 



B.2.2. Linearization 

Here we approximate Eqs. (jT4j» to first order for small pertur- 
bations around the hydrostatic equilibrium states. We substitute 
Eqs. (jT3J in Eqs. ( |T4~] > and expand to first order in quantities with 
subscript 1. By noting that M(ro, t) = M(ro, t = to), we obtain: 



dpi p Q d 2 



dt 



dr 



(M = , 



dui kT dpo/dr ( n dpjdr TA GM(r ) 
12— + - — 7^— + — 1 + 2 — r\ 



dt 



m p 



r dpo/dr Q T 



From the usual rules of derivation, we have: 
dr(r , t) dr (r, t)/dt 



dt 



dro(r, t)/dr 



(B.ll) 



(B.12) 



By applying Eqs. ( |B.10| i and ( |B.8| >, Eq. ( |B.12| > can be written as: 

5-="i. (B.13) 
dt 

Now we wish to obtain one equation involving only the un- 
known pi, by combining the three Eqs. |BTTT]i and ( |B.13|l. We 
assume the modal dependence ( |16) , We eliminate r t from ( |B. 13| 
and the second of ( |B.ll| i. We then eliminate u\ from the re- 
sulting equation and the first of ( |B.l \\ , to obtain the follow- 
ing equation {we also used hydrostatic equilibrium to replace 
[(dp /dr )/p ](kT /m) = -[GM(r )/r 2 Q ]}: 



rl dr 



I GM(r ) 



= 0. 



(B.14) 



By referring to the dimensionless radius £o = fo/A, we obtain 
Eq. ([17]). 

6.3. Temperature expression for the constant {E, V} case 

Here we show the steps leading from Eq. ([26]) to Eq. p7] ). 

From the Poisson and Emden equations, the dimensionless 
potential if/ is related to the gravitational potential Oo of the un- 
perturbed density distribution in the following way: 



W) = [fcoW - *o(0)] KkT /m) 



(B.15) 



From Eq. ( |B.15| l, by recalling that we consider only perturba- 
tions that do not change the total mass ( Jp!(r)d 3 r = 0), we 
obtain: 

f Pl (r)^(r)d 3 r 
Ti = 7- T . 



(B.16) 



From Eq. ( |B.3| ), by setting / = pou\ and £ = r/A, we obtain 
Eq. ([27]). 

Appendix C: Detailed angular momentum 
conservation and collapse 

Here, for completeness, we show in detail that the angular mo- 
mentum barrier prevents a spherically symmetric collisionless 
system from collapsing. We consider only perturbations that do 
not break the assumed spherical symmetry of the system. 

Consider a collisionless system of particles of individual 
mass m and total mass M. From the assumption of spherical 
symmetry, each particle is confined to a plane and we can write 
the single-particle Lagrangian as: 



m(r 2 + r 2 ^) - V(r, t) 



(CI) 



where V(r, f) is a time-dependent potential, t is time, r is the 
distance from the center and is the angular coordinate. The 
Lagrangian dC 1| > conserves the angular momentum of the parti- 
cle, that is, r 8 - C, where C is a constant. Then the equation of 
the motion is: 



C 2 1 dV(r,t) C 2 GM(r,t) 



(C2) 



where M(r, t) is the total mass contained in the sphere of ra- 
dius r. Equation \C2\ is the equation of the motion of a par- 
ticle moving in one dimension and subject to the force F r = 
mC 2 /r 3 - GmM(r, t)/r 2 . The following inequality holds: 



< M(r, t) < M . 



(C3) 
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By multiplying by r and integrating both sides of Eq. \C2) , we 
obtain: 



rHt) 



nt ) { c 2 

2 + 2 



1 



1 



r z (t ) r\t) 



rr(t) 
Jr(t ) 



GM(s, t(s)) 



ds. 

(C4) 

Since r 2 (t)/2 is positive, the right-hand side of Eq. ( |C.5| l must 
be positive. By taking r(t) < r(fo) (we are not interested in the 
case in which r(t) is greater than the initial radius) and using 
Eqs. (|C.5|) and (|C 3|>, we find: 



r 2 (t) 



- f lM | C 

2 2 

< , C2 

"2 2 

2 2 



1 1 

1 1 



r 2 (f ) r 2 (t) 
1 1 



X 

f 



GM(s,t(s)) 



ds 



'•('o) Q M 



-V) 
+ GM 



ds 

1 1 

r(t) r(t ) 



(C5) 



For given values of r(?o) and r(fo), the quantity appearing in the 
last line of Eq. |C3) tends to -oo as r(t) — > 0. Hence, for given 
initial conditions the particle cannot reach arbitrarily small val- 
ues of r(t). 



Appendix D: Density and velocity profiles of the 
linear modes 

In this appendix we show density and velocity profiles of the nor- 
mal modes for the linear stability analysis presented in Sect. |4j 
Pi/po and u\ are meant to be in arbitrary scales. 
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DA. Constant {T,V} profiles 

L = -0.02 L = -0.02 




Fig. D.l. Relative density perturbation profiles pi(£)/po(£) and velocity profiles of the normal modes for the constant {T, V} case, 
obtained by solving Eq. ( fT8| ). The profiles should be truncated at a value ^ = S where the velocity profile vanishes, to satisfy 
boundary conditions ( fTTj ). The vertical dotted line indicates where the system should be truncated to obtain the mode of lowest L 
for fixed H: for L = -0.02 and L = only this mode is entirely displayed, while for L = 0.02 two modes are displayed, depending 
on which zero of the velocity profile is chosen. In the case L = 0, the total density p(f) = po + Pi(0 at the point £ = 4.07 remains 
unchanged, that is, unperturbed; this is one of the relevant points listed by |Lynden-Bell & Wood ( 1968| l. 
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D.2. Constant {E, V] profiles 

L = -0.02 L = -0.01 




-0.2- : -0.2 



20 40 60 80 20 40 60 80 




-0.2- ! -0.2 



10 20 30 40 10 20 30 40 




10 20 30 40 10 20 30 40 



Fig. D.2. Relative density perturbation profiles pi(£)/po(£) of the normal modes for the constant {E, V} case, obtained by solving 
Eq. ( f28| . The modes should be truncated at a value £ = S where the corresponding velocity profile shown in Fig. |D.3| vanishes, as 
marked by the vertical dotted lines, in order to satisfy the boundary conditio ns (fTTj) . Zeros are displayed. Only modes of lowest L at 
given E are shown. Note that the core-halo structure described in Subsection |4. 1.2| disappears between L = 0.021 and L = 0.022. 
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20 



Fig. D.3. Velocity profiles of the normal modes for the constant [E, V} case, obtained by solving Eq. ( |28| ). The modes should be 
truncated at a value £ = H corresponding to the vertical dotted lines, in order to satisfy the boundary conditions ( fTT) . Other zeros are 
displayed. Only modes of lowest L for fixed S are shown. Note that the core-halo structure described in Subsection 4.1.2 disappears 
when the velocity has no internal zeros, between L - 0.021 and L - 0.022. 
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D.3. constant {T,P} profiles 

L = -0.03 L = -0.02 




Fig. D.4. Relative density perturbation profiles pi(£o)/Po(£o) of the normal modes for the constant {T, P] case, obtained by solving 
Eq. ( |38| >, calculated and displayed here in the Lagrangian representation. The modes should be truncated at a value £o = 3 where 
the density profile vanishes, in order to satisfy the boundary conditions ( [39] ). The first zero, which represents the mode of minimum 
L at given H, is indicated by the vertical dotted line. A mode of higher L for fixed S is shown in the L = 0.03 case; for other cases 
higher modes can be identified in a similar way. 
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Appendix E: Equations for the two-component case 

In this appendix we summarize the equations of the linear anal- 
ysis for the two-component case and show some examples of 
the density profiles associated with the modes that characterize 
the onset of the instability. We denote by subscripts A and B the 
lighter and the heavier component, respectively. 

The unperturbed states are the two-component self- 
gravitating truncated isothermal spheres considered byjTafFeFaT 
([T975)l; pghtrnanl ([I977l>; |Yoshizawa et al.| ( [T978) ; |de Vega & 
Siebert ( 2002| l; Sopik et al. ( 2005 1. The density profiles can be 



written as: 



Pao(0 = Pao(OK 
Pm(£) = Pm(0)e 



-U+W 



Pao(0 
.Pbo© 



if^<E 



iff >S, 



(E.l) 



where pao and p B o are respectively the density profiles of the 
lighter and heavier component, f = rjXq. is the dimensionless 
radial coordinate, where X% = [kT(l/m A + l/m B )/4-7rGpo(0)] l/2 
and we denote by po(f ) = Pao(0 + Pbo(^) the total unperturbed 
density; S is the value of f at the truncation radius, B = m B /m A 
is the ratio of the single-particle masses, ifr is the solution of the 
following generalization of the Emden equation |5j: 



df 



(*V) = f 



1 + a 



1 + 



m = 0/(0) = 0, 



(E.2) 



(E.3) 



where a = Pao(0)/Pbo(0) is the ratio of the unperturbed central 
densities. The symbol ' denotes derivative with respect to the 
argument f . 

The linearized hydrodynamical equations, governing the 
evolution of the two-component fluid system for small devia- 
tions from the unperturbed states described above, are obtained 
by generalizing in a straightforward manner the steps leading 
from Eqs. Q to Eq. ( fT8} , The result, which generalizes Eq. ( fTS) , 
is the following system of equations that governs the evolution 
of radial perturbations: 



while the two conditions at the truncation radius satisfy the re- 
quirement that the radial velocities must vanish at the edge. 

The system ( |E.4| i for L — is equivalent to the system that 
can be obtained by generalizing in a straightforward manner the 
thermodynamical analysis of Chavanis ( 2002| l. The latter analy- 
sis can be used to find the points for the onset of instability. This 
proves that the onset of instability occurs at the same values of 
5 in the dynamical and in the thermodynamical approach. 

In Fig. E.l we show the density profiles for the marginally 
stable modes (L - 0) in three different situations, that is, with 
B — 3 and three different values of Mb/Ma- The density per- 
turbation of the heavier component is greater than the density 
perturbation of the lighter component even for small values of 
Mb I Ma, indicating that the heavier component is the more im- 
portant driver of the instability. 



L/a 



1 + 



P 



V\ffA+f A 



2 f 2 f 
~^JA + -^JA 



1 + 



1 + 



-e-^(fA+fs) 



Lf B = 



- (1 +B) ij/ + - f' B ' - 2 -f B + ^f B 

-—*- (1+m (fA+f B ). 



l+B 



1 + a 



(E.4) 

Here L = w 2 /4^Gpo(0) represents the dimensionless (squared) 
eigenfrequency, / A (f) = p A0 (f)MAi(f) and / fi (f) s p B0 (f)Wi?i(f), 
where uai and ub\ are the radial velocity perturbations of the two 
components. The boundary conditions are: 



/a(0) = /b(0) = 
f A (E) = f B (S) = . 



(E.5) 



Similarly to the one-component case, the two conditions at the 
center follow from requiring regularity and spherical symmetry, 
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{L=0, a=10, j8=3, M B /M A =1.78) {L=0, a=l, /J=3, M B /M A =0.066) 





{L=0, ff=0.1, j8=3, M a /M /4 =0.0008) 



Fig. E.l. Relative density perturbation profiles P\a(^)/Po(£) (lighter component, dotted line) and P\b(£)/Po(^) (heavier component, 
solid line) of the normal modes for the two-component constant {T, V} case, obtained by solving Eq. ( |E.4| >. The plots show marginally 
stable modes (L = 0) of minimum L at given S. They represent the density profiles that characterize the onset of the instability for 
fixed value of /3 — m B /m A = 3 at different values of the total mass ratio M B /M A . Even for small values of M B /M A , the density 
perturbation of the heavier component dominates, thus suggesting that the heavier component is the more important driver of the 
instability. 
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